#include<iostream>
#include<cmath>

#include "global.hpp"
#include "solver.hpp"

using namespace std;
//typedef enum {
//  CGNS_ENUMV( BCTypeNull ) =CG_Null, used as nonBC
//  CGNS_ENUMV( BCTypeUserDefined ) =CG_UserDefined, used as BCOverset
//  CGNS_ENUMV( BCAxisymmetricWedge ) =2,
//  CGNS_ENUMV( BCDegenerateLine ) =3,
//  CGNS_ENUMV( BCDegeneratePoint ) =4,
//  CGNS_ENUMV( BCDirichlet ) =5,
//  CGNS_ENUMV( BCExtrapolate ) =6, used as PML
//  CGNS_ENUMV( BCFarfield ) =7, 
//  CGNS_ENUMV( BCGeneral ) =8,
//  CGNS_ENUMV( BCInflow ) =9, 
//  CGNS_ENUMV( BCInflowSubsonic ) =10,
//  CGNS_ENUMV( BCInflowSupersonic ) =11,
//  CGNS_ENUMV( BCNeumann ) =12,
//  CGNS_ENUMV( BCOutflow ) =13, 
//  CGNS_ENUMV( BCOutflowSubsonic ) =14,
//  CGNS_ENUMV( BCOutflowSupersonic ) =15,
//  CGNS_ENUMV( BCSymmetryPlane ) =16,
//  CGNS_ENUMV( BCSymmetryPolar ) =17,
//  CGNS_ENUMV( BCTunnelInflow ) =18,
//  CGNS_ENUMV( BCTunnelOutflow ) =19,
//  CGNS_ENUMV( BCWall ) =20,
//  CGNS_ENUMV( BCWallInviscid ) =21,
//  CGNS_ENUMV( BCWallViscous ) =22,
//  CGNS_ENUMV( BCWallViscousHeatFlux ) =23,
//  CGNS_ENUMV( BCWallViscousIsothermal ) =24,
//  CGNS_ENUMV( FamilySpecified ) =25, used as BCImpedance
//} CGNS_ENUMT( BCType_t );

void Solver::BCFarfield_Radiation(double tRK, Mesh &msh, Node *node,Solution *sol,Flux *flux, Meanflow *mf){
    int iZone,iLZone,iLNode,iDim,iDimC,iFace,iEq,dimBC;
    int ijk[NDIM],iSt[NDIM],iEd[NDIM];
    int dN[NDIM],dI,iT,iS;
    double dQC[NDIM];
    double cSquare,Vtheta,xyz[NDIM],rtp[NDIM];
    double Q[NEQ],dQ[NDIM],MQ[NEQ],dQr;
    double MU[NDIM],uVecR[NDIM],uVecTheta[NDIM],uVecPhi[NDIM];
    //Calculate the flux of the points which are in the wall regions
    for (iLZone=0;iLZone<msh.nLZone;++iLZone){
        iZone=msh.lZoneList[iLZone];
        dN[0]=1;
        dN[1]=msh.size[iZone][0];
        #ifndef _2D
        dN[2]=msh.size[iZone][0]*msh.size[iZone][1];
        #endif
        for (iFace=0;iFace<NFACE;++iFace){
            if (msh.bcTyp[iZone][iFace]==7){
                for (iDim=0;iDim<NDIM;++iDim){
                    iSt[iDim]=msh.bcRng[iZone][iFace][iDim][0];
                    iEd[iDim]=msh.bcRng[iZone][iFace][iDim][1];
                }
                #ifndef _2D     
                for (ijk[2]=iSt[2];ijk[2]<=iEd[2];++ijk[2])
                #endif
                for (ijk[1]=iSt[1];ijk[1]<=iEd[1];++ijk[1])
                for (ijk[0]=iSt[0];ijk[0]<=iEd[0];++ijk[0]){
                    iLNode=msh.caliLNode(iZone,ijk);
                    for (iDim=0;iDim<NDIM;++iDim){
                        xyz[iDim]=(node+iLNode)->xyz[iDim];
                        MU[iDim]=(mf+iLNode)->MQ[1+iDim];
                    }
                    cSquare=(mf+iLNode)->MC*(mf+iLNode)->MC;
                    #ifndef _2D
                    math.CartesianToSpherical(xyz,rtp,uVecR,uVecTheta,uVecPhi);
                    Vtheta=math.vectorDotProduct(NDIM,MU,uVecR)+sqrt(cSquare-pow(math.vectorDotProduct(NDIM,MU,uVecTheta),2)-pow(math.vectorDotProduct(NDIM,MU,uVecPhi),2));
                    #else
                    math.CartisianToPolar(xyz,rtp,uVecR,uVecTheta);
                    Vtheta=math.vectorDotProduct(NDIM,MU,uVecR)+sqrt(cSquare-pow(math.vectorDotProduct(NDIM,MU,uVecTheta),2));
                    //cout<<ijk[0]<<", "<<ijk[1]<<": "<<rtp[0]<<" "<<rtp[1]<<" "<<Vtheta<<endl;
                    #endif
                    for (iEq=0;iEq<NEQ;++iEq){
                        //Calcluate   derivative in computational domain
                        for (iDimC=0;iDimC<NDIM;++iDimC){
                            dQC[iDimC]=0;
                            iT=(node+iLNode)->locTyp[iDimC];
                            for (iS=0;iS<STENCIL;iS++){
                                dQC[iDimC]+=s.d[iT][STENCIL-iT-1][iS]*(sol+iLNode-(iT-iS)*dN[iDimC])->Q[iEq];
                            }
                        }
                        for (iDim=0;iDim<NDIM;++iDim){
                            (flux+iLNode)->dQ[iEq][iDim]=0;
                            for (iDimC=0;iDimC<NDIM;++iDimC){
                                (flux+iLNode)->dQ[iEq][iDim]+=dQC[iDimC]*(node+iLNode)->ja[iDimC][iDim];
                            }
                            dQ[iDim]=(flux+iLNode)->dQ[iEq][iDim];
                        }
                        dQr=math.vectorDotProduct(NDIM,dQ,uVecR);
                        (flux+iLNode)->FRK[iEq]=-Vtheta*(dQr+(sol+iLNode)->Q[iEq]/2.0/rtp[0]);
                    }
                }
            }
        }
    }      
}

void Solver::BCFarfield_Outflow(double tRK, Mesh &msh, Node *node,Solution *sol,Flux *flux, Meanflow *mf){

}

void Solver::BCFarfield_PML(double tRK, Mesh &msh, Node *node,Solution *sol,Flux *flux){
    
}

